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Abstract 

We consider a plane channel flow of an electrically conducting fluid which is driven by a mean pressure 
gradient in the presence of an applied magnetic field that is streamwise periodic with zero mean. Magnetic 
flux expulsion and the associated bifurcation in such a configuration is explored using direct numerical 
simulations (DNS). The structure of the flow and magnetic fields in the Hartmann regime (where the 
dominant balance is through Lorentz forces) and the Poiseuille regime (where viscous effects play a significant 
role) are studied and detailed comparisons to the existing one-dimensional model of Kamkar and Moffatt 
(J. Fluid. Mech., Vol.90, pp 107-122, 1982) are drawn to evaluate the validity of the model. Comparisons 
show good agreement of the model with DNS in the Hartmann regime, but significant diferences arising 
in the Poiseuille regime when non-linear effects become important. The effects of various parameters like 
the magnetic Reynolds number, imposed field wavenumber etc. on the bifurcation of the flow are studied. 
Magnetic field line reconnections occuring during the dynamic runaway reveal a specific two-step pattern 
that leads to the gradual expulsion of flux in the core region. 
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1. Introduction 


Interaction of electrically conducting flows with magnetic fields which forms the basis of the subject 
of magnetohydrodynamics (MHD), occur commonly in nature and are frequently utilized in industrial pro¬ 
cesses. Naturally they occur in astrophysical phenomena like the formation of stars and galaxies, phenomena 
in the Sun like sunspots and solar flares, and in the generation of terrestrial magnetic field in the core of 
the Earth by the dynamo action [l| . In industrial processes like continuous casting of steel and aluminium, 
material processing and plasma confinement in nuclear fusion applications, magnetic fields are applied ad¬ 
vantageously for flow manipulation and control M. In all MHD flows, the magnetic field affects the flow 
through Lorentz forces, but the back reaction of the flow on the magnetic field (that leads to the bending 
of field lines) depends on a parameter called the magnetic Reynolds number defined as. 
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where U and L are the charateristic velocity and length scales in the flow and A is the magnetic diffusivity 
of the fluid given by A = (/iocr)~^, l^o and a being the magnetic permeability of free space and the electrical 
conductivity of the fluid respectively. Significant bending of magnetic field lines occur only when Rm ^ 1 
and larger. Flows considered in this paper fall into this category. 

An interesting feature of MHD flows at high magnetic Reynolds numbers is the expulsion of magnetic flux 
that typically occurs under the imposition of an electrically conducting fluid flow with closed streamlines. 
This follows from an analogy of the well known Prandtl-Batchelor theorem Q in classical hydrodynamics. 
The kinematic problem of magnetic flux expulsion under rotation has been extensively studied during the 
sixties (see for e.g. [1-0) in the context of astrophysics. That flux expulsion also persists in the dynamic 
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regime was pointed out by Galloway et al. [9[ and was followed by further analysis of the dynamic effects of 
flux ropes in Rayleigh-Benard magnetoconvection IS- 

An important aspect of the dynamic behavior associated with flux expulsion is the ‘runaway’ effect, 
which can be explained as follows. When magnetic field lines start to get expelled in a region of the flow, 
there is a decrease in the Lorentz forces that opposes the mean flow leading to an acceleration of the flow in 
that region. This in turn leads to further expulsion of magnetic flux and subsequently results in a cascading 
effect of flow acceleration and flux expulsion, wherein dissipative forces like viscosity ultimately balances 
the driving force, leading to a steady state. This effect can play a significant role in in the performance of 
electromagnetic pumps that are used to pump liquid metal. Early analytical studies of this phenomenon by 


Gimblett et al. [ll| using rotating cylindrical and spherical solid bodies under an applied normal magnetic 


field showed an associated Thom cusp catastrophe and hysteresis effect. Such behaviour was seen further in 
the fluid context by Moffatt [l2| . 

However, flux expulsion can also happen in flow configurations without closed streamlines, if the imposed 
magnetic field is non-uniform and periodic in the mean flow direction. A particularly interesting configuration 
is that of a plane channel flow driven by a mean messure gradient with an imposed sinusoidal magnetic 
field that was analysed by Kamkar and Moffatt [l3| which will be denoted as KM82 hereafter. In their 
study, the interaction of the flow and magnetic fields was described by simplified one-dimensional model 
equations. Steady state solutions were obtained from which two different flow regimes were identified, 
namely the Hartmann and Poiseuille regimes and the location of the bifurcation leading to the transition 
between these two regimes was computed. However, various simplifications were assumed in that study. 
For example, the non-linear terms (and hence the Reynolds stress terms) in the Navier-Stokes equation 
and the variations along the streamwise direction were neglected. Although it enables one to obtain quick 
solutions, the approximate model can lead to significant loss of accuracy and underprediction/overprediction 
of the jump that occurs during the bifurcation. The focus of our work is to perform 2D direct numerical 
simulations of the problem similar to KM82 with a twofold purpose. On one hand, it helps validate the ID 
model predictions at the steady state and quantify the differences arising out of the simplifications of the 
model. On the other hand, the presence of non-linearities can result in time-dependent solutions for the 
flow and magnetic fields in both the regimes in the final state. 

The paper is organized as follows. Section [2] describes the problem setup and the full governing equations 
along with a brief overview of the KM82 model. This is followed by the details of the numerical procedure 
in section [3l In section HI numerical results of the DNS are presented and compared to the predictions of 
the model, including the effect of various parameters on the charateristics of the bifurcation, followed by 
conclusions in section [5l 


2. Problem setup and governing equations 

2.1. Problem formulation and full governing equations 

We consider the two-dimensional incompressible flow of an electrically conducting fluid (e.g. a liquid 
metal) driven by a mean pressure gradient in a straight rectangular channel. Periodicity is assumed along 
the streamwise direction x and the wall normal direction is denoted by z. A magnetic field Bo{x,z) = 
B{x,z,t = 0) (generated by electric current or magnet sources outside the channel) with a prescribed 
wall normal component Bqz = cos(kx) with a wavenumber k, is imposed on the flow. Such a periodic 
magnetic field can be ideally produced for e.g. by magnets distributed on the channel walls with alternating 
north and south poles in the streamwise direction (see 0 ). The action of the flow on the magnetic field 
generates plane-normal electric current densities J = (0, J(x, z), 0) in the flow which leads to the generation 
of a secondary magnetic field and Lorentz forces that affect the flow. The imposed magnetic field which 
is divergence-free can be expressed as Bq = V x Aq, where Aq = (0, Ao(x, 2;), 0) is the magnetic vector 
potential. Ghoosing the scales of half channel height L for the length, the maximum value of the imposed 
magnetic field Bq for the magnetic field, BqL for the vector potential and applying the curl-free condition 
on Bo, we get in the non-dimensional form, 

9^ao d‘^ao 
dx‘^ ^ dz^^ 
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with the boundary conditions 
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= cos{nx) 


on z = ±1 ; ao (0, z) = ao {lx, z) 


( 3 ) 


where the symbol n = kL represents the normalized wavenumber, is the non-dimensional streamwise 
length of the channel and ao is the normalized vector potential. Solution of equation (|2]) using seperation 
of variables yields 

1 sin(Kx)cosh(Kz) 

ao =- w ^ • 4 

K cosn{n) 

This leads to the form of the non-dimensional initial (or imposed) magnetic field 6o as 


sin{Kx)sinh{Kz) . cos{Kx)cosh{Kz) 

Oo = - ^ H-- k. 


cosh{i<i) 
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( 5 ) 


the field lines of which are shown in Fig. [H Here i and k refers to the unit vectors in the streamwise {x) 
and wall-normal {z) directions respectively. 



X 


Fig. 1. Field lines of the imposed magnetic field 6o as given by equation Two magnetic X-points at 
{x = 'K{2^ z = 0) and {x = Z'k j2^ z — 0) can be observed at the centerline. 

The physics of the problem is governed by the Navier-Stokes equation for the momentum balance in¬ 
cluding the additional source term representing the Lorent force (body force) produced by the induced 
electric currents and the induction equation for magnetic field transport along with the constraints of mass 
conservation (continuity) and solenoidality of the magnetic field. The pressure gradient is decomposed as 
VPt = —pGi + VP, where VPt represents the cumulative pressure gradient and —pG the constant mean 
pressure gradient that is applied in the streamwise direction. Non-dimensionalizing using the scales 
X/L‘^kG^ pGL and crXBolL‘^k for velocity, time, pressure and the current densities respectively, and denoting 
all non-dimensional variables by small letters, the system of governing equations can be written as 
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^ + ^(^.V)^=l_Vp + «V»„ + iox6), 

V • V = 0, 

V-6 = 0, 


j = K (V X 6), 

u = w = 0, bz = cos{kx) on ^ = ±1 ; v,b periodic in x-direction 
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with the no-slip and no penetration boundary conditions for the fluid velocity on the walls and periodicity 
assumed in the streamwise direction. Furthermore, the wall normal component of the magnetic field hz 
on the walls remain unchanged (equal to the imposed magnetic field, h^z) and the streamwise component 
follows from the divergence-free condition The parameters involved in the problem are 


_ v\ L^kG _ pL^kG 


n = kL 


( 12 ) 


where u and p represent the kinematic viscosity and the mass density of the fluid respectively. The parameter 
[3 represents the magnetic Reynolds number and the parameters e, Q can be regarded as the inverse of the 
hydrodynamic Reynolds number and the Stuart numbers respectively. All the fluid properties are assumed to 
be constant. The coupled evolution of the velocity and magnetic fields is computed by solving the governing 
equations numerically, a brief summary of which is presented in the next section. 


2.2. One-dimensional approximate model of Kamkar and Moffatt (KM82) 

In view of the fact that the results of our DNS are compared with the ID approximate model of KM82, 
we present a brief overview of their model here. The following assumptions have been made in the model 

• The secondary magnetic field consists of only a single mode (wavenumber) along the streamwise 
direction, which is taken to be the wavenumber {k) of the applied magnetic field. To this effect, the 
vector potential A of the magetic field h is expanded as 

A(x, z, t) = (13) 

where 5R represents the real part and f{zA) is the dimensionless profile function. 

• The variation of dynamics in the streamwise direction is neglected and hence the mean (x-averaged) 
governing equations are considered. 

• The velocity fluctuations v' are small compared to the mean (x-averaged) streamwise velocity I/, 
jv'l ^ U and hence the Reynolds stress terms in the momentum equation and the fluctuating parts 
of advection terms in the A-transport are neglected. This assumption is supposed to be valid when 

The assumptions stated above lead to the following mean governing equations for U{zA) and f{zA) in the 
non-dimensional form. 
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f/ = 0, / = 1 on 2 = ±1 ; 


4 


(15) 

(16) 








which correspond to equations (2.33) and (2.38) in KM82, where /* is the complex conjugate of /. For later 
comparisons with the DNS results, we solve the above model equations ([T4|) to ([T 6 ]) by a finite-difference 
method on a uniform grid. 


3. Numerical Procedure 


The numerical solution of the system (| 6 |) to (pT]) is obtained using a second order finite difference method. 
The domain is discretized into a rectangular Cartesian grid with uniform grid spacing along the streamwise 
direction and a non-uniform stretched grid in the wall normal direction in order to resolve the thin Hartmann 
boundary layers near the walls, that are typical of wall bounded MHD flows. The stretched grid is obtained 
by a coordinate transformation from a uniform grid coordinate 77 to the non-uniform grid coordinate z 
through 


tanh(S'77) 

tanh(5') 




(17) 


where S represents the stretching factor. A typical grid used in our studies is shown in Fig. [2l 



Fig. 2 . Non-uniform structured grid in the x^-plane with 129 x 129 grid points and S = 2 . 2 . 


The momentum equation (| 6 ]) is integrated by the standard projection scheme where the non-linear, 
viscous and Lorentz force terms are treated explicitly by the Adams-Bashforth method to first compute an 
intermediate velocity field v* which is then projected onto a solenoidal velocity field at the new time 

level by a correction obtained from the solution of a Poisson equation for the pressure. A second order 
backward finite difference scheme is used for time discretization using the levels n — 1 , n, n 1 where n 
is the current time level. Details of the procedure to compute the velocity field can be found in Krasnov 
et al. HI. For the magnetic field, the normal component of the magnetic field bz is computed by using 
a semi-implicit procedure, wherein only the diffusive term is treated implicitly and the advective and field 
stretching terms are treated explicitly. This requires again the solution of a Poisson equation for bz. The 
FORTRAN software package FISHPACK [l5| has been used to solve the Poisson equations. Subsequently, 
the streamwise component bx is reconstructed from bz using equation m, in order to satisfy the solenoidality 
of the magnetic field. Alternatively, it is possible to solve for the vector potential A (since the problem is 
2 D, A has only one component) and recover the magnetic field components bx and bz from A. The magnetic 
field so obtained is used to compute the current density j which is required to evaluate the Lorentz force 
term in the Navier-Stokes equation. OpenMP parallelization has been used in performing the computations 
for the results presented in this paper. 
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4. Results and comparison 

Starting with either an initial laminar velocity profile (with no streamwise variation) or fluid at rest 
{v = 0) and the imposed magnetic field 6o, the governing equations are numerically integrated in time to 
obtain the final equilibrium states. All the computations have been performed for a streamwise domain 
length of one period, lx = 27r/K on a 129 x 129 grid. A grid sensitivity study was performed which indicated 
that further increase in grid resolution does not improve the solution accuracy significantly, within the 
parameter space studied in this paper. Depending on the velocity profiles of the final states, two regimes of 
flows are defined (as in KM82), namely the Hartmann regime and the Poiseuille regime. 




Fig. 3. Steady state streamwise velocity profiles in the Hartmann regime at (a) x = 7r/2 and (b) x = n. 
Parameters are (3 = k, = 1 and e = 5 x 10“^. 

The Hartmann regime is characterized by very steep velocity gradients in the boundary layers as com¬ 
pared to the core region. Pressure gradient in the core is dominantly balanced by Lorentz forces whereas in 
the boundary layers it is a combination of viscous and Lorentz forces that balances the pressure gradient. 
Flows at relatively small Q or high interaction parameter belong to this regime. In contrast, the Poiseuille 
regime typically demonstrates Toiseuille-like’ (parabolic) axial velocity profiles and is dominated by viscos¬ 
ity in the core region. Flows at relatively higher Q belong to the Poiseuille regime. The final equilibrium 
state of the Hartmann regime is observed to be steady in time unlike the Poiseuille regime where significant 
velocity fluctuations persist in final state. Transition between the two regimes occurs over a narrow band of 
Q through a bifurcation that is a manifestation of the runaway effect. We now provide a very brief account 
of the nature of the solutions in these two regimes. 

4 . 1 . Hartmann regime 

Typical axial velocity profiles in the Hartmann regime are shown in Fig. [3] at two different axial locations 
of the channel, x = ^ and x = tt and at various values of the parameter Q. These two locations correspond 
to the streamwise extreme values of the imposed magnetic field bo. It can be observed from Fig. [3] that 
higher axial velocities (or flow acceleration from the initial state) are observed with increase in Q and the 
profiles in the boundary layers at x = tt look more ‘Hartmann-like’ than at x = ^ due to the pronounced 
effect of the wall normal component boz. The distribution of axial velocity Vx in the domain as shown 
in Fig. |4(a)| (along with the velocity streamlines) clearly indicates the laminar nature of the flow in the 
Hartmann regime. Advection of the magnetic field can be observed from the corresponding field lines shown 
in Fig. |4(b)| that indicate only a slight bending of the field lines. As is clear from the velocity field, no 
significant events of severing or reconnection of field lines are observed in this regime. 
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(b) 

Fig. 4. Typical streamwise velocity and magnetic field configuration in the Hartmann regime, (a) Contours 
of streamwise velocity and (b) magnetic field lines in the steady state at Q = 0.3, for p = 1, n = 1 and 
6 = 5x10“^. We observe how the magnetic X-points are slightly shifted by the flow. 
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(c) 


Fig. 5. Contours for the evolution of velocity field to the Poiseuille-like state as the flow accelerates from 
rest, (a) t = 0, (b) t = 16.1, (c) t = 60.1. Parameters are Q = 0.5, {3 = n = 1 and e = 5 x 10“^. Left 
column: streamwise velocity Vx- Right column: wall-normal velocity Vz- 


4 . 2 . Poiseuille regime 

The Poiseuille regime is a result of the runaway effect, where there is a significant acceleration of the 
flow due to considerable bending and severing of field lines and subsequent expulsion of magnetic flux. The 
streamwise velocity profiles show significant gradients in the wall normal direction and hence look ‘Poiseuille- 
like’. It is observed that the flow exhibits strongly unsteady behaviour even in the final (steady on average) 
state. Of particular interest is the initial phase of the transient flow that ensues when the flow starts to 
accelerate from rest, due to the applied mean pressure gradient. The evolution of the velocity field leading 
to almost complete expulsion of magnetic flux in the core {z = 0) in such a case is shown in Fig. [5] through 
snapshots of velocity component contours. It can be seen that the gradual acceleration of the streamwise 
velocity is accompanied by relatively small wall normal velocity component on both side sides of the core, 
in a staggered arrangement. At the same time, the normal component of the magnetic field (which leads to 
streamwise Lorentz force) in the core is gradually destroyed, due to the bending of the vertical field lines 
(due to advection) as well as the reconnection of the field lines, depending on the streamwise location in 
the channel. Reconnection here refers to the rearrangement of magnetic field line topology. This can be 
observed in Fig. [6l where snapshots of the configuration of the magnetic field lines show the eventual decay 
of magnetic flux in the core and significant shifting (advection) of the X-points, due to reconnection events 
(the mean flow is from left to right). A specific common pattern in the reconnetion of the magnetic field 
lines was observed during the process of dynamic runaway. A field line in the region with a strong negative 
wall normal component hz (those that are approximately located at 2.5 < x < tt in Fig. |6(a)[ ), undergo a 
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(c) 

Fig. 6. Advection and expulsion of magnetic flux as ^he flow starts to accelerate from rest, (a) t = 0, (b) 
t = 16.1, (c) t = 60.1. Parameters are Q = 0.5, /3 = 1, /i: = 1 and e = 5 x 10“^. Contours coloured by hz^ 



































































Fig. 7. Severing and reconnection of magnetic field lines shown in a small section of the channel when the 
flow accelerates from rest, (a) t = 0, (b) t = 2.0, (ff) ^ = 2.1, (d) t = 2.9, (e) t = 3.0, (f) t = 15.4, (g) 
t = 15.5, (h) t = 37.0. Parameters are Q = 0.5, ^ n = 1 and e = 5 x 10“^. Following the line marked 
AA’ shows a characteristic reconnection and stretching pattern that it undergoes leading to the expulsion 
of magnetic flux in the core. 
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Fig. 8. Evolution of the wall normal component of the magnetic field (bz) along the channel centerline z = 0 
when the flow accelerates from rest (obtained from DNS). Parameters are /3 = 1, Q = 0.5 and e = 5 x 10“^. 
Inset shows corresponding time decay of the mean magnetic energy on the centerline < > at 2 ; = 0. We 

mark the magnetic X-points and their steady shift further downstream by reconnection events. 


two-fold reconnection process and transform into a field line with positive (except at the core, where it 
is zero). A typical example of this is shown in more detail in Fig. [71 where the field line (marked at the 
ends by A and A’) corresponding to the imposed magnetic field initially develops a sharp ‘pinch’ at 2 ; = 0 
(Fig. |7(b)| ) before a reconnection event leading to both ends of the line attached to the top wall (Fig. |7(c) ). 
After some further stretching, another reconnection event occurs as seen from Fig. |7(d)] to Fig. |7(e)' leading 
to a reversal of the direction of the magnetic field as compared to the initial state. Subsequently the field 
line is stretched significantly in the flow direction as shown from Fig. |7(d)| through Fig. |7(h)| leading to the 
final topology of the line that also results in 6;^ = 0 at 2 ; = 0. Interestingly, these series of events is seen to 
occur to every flux line (in the region considered, 2.5 < x < tt) in a sequential manner from left to right. 
This is clearly seen from the pinching and reconnections occuring to the field line next to AA’ seen from 
Fig. |7(f) I through Fig. |7(h)[ 

The decay of the magnetic flux in the core during this period (from t = 0 to t = 370) is shown in Fig. (HI 
where the evolution of bz(x) along the channel centerline is plotted. Advection of X-points in the streamwise 
direction and the decay of the amplitude of b^ can be clearly observed. This is accompanied by the temporal 
decay of the mean magnetic energy at the channel centerline which defined as 


< Ef, >= 




(18) 
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Fig. 9. Snapshots of streamlines of the velocity field in the final (steady on average) state showing the 
transport of vortices near the wall. Coloured by contours of (a) t = 0.3, (b) t = 1.1, (c) t = 1.9. 
Parameters are Q = 0.5, /3 = 1, k. = 1 and e = 5 x i(j^. 














































































































(b) 


Fig. 10. (a) Steady state mean streamwise velocity profiles in the Hartmann regime, (b) Mean streamwise 
velocity profiles in the Poiseuille regime (streamwise and time averaged). Parameters are /3 = 1, /^ = 1 and 
e = 5 X 10“^ 


This decay contributes to the growth of kinetic energy and hence the runaway process. 

As mentioned previously, the final state occuring in the Poiseuille regime shows a strongly unsteady 
behaviour, with secondary flow structures which become more frequent at lower viscosities. A typical 
example is shown in Fig.[9l where two large vortex structures are observed, one on either wall and seperated 
in the streamwise direction. These vortices (or recirculation zones) are advected along the mean flow 
(Figs. |9(a)| through 9(c)[ ), leading to an almost time-periodic behaviour of the flow. Absence of chaotic 
states might be attributed to low Reynolds numbers and short domain length in the problem. 


4 . 3 . Comparison with the predictions of KM82 

We now turn to the comparison of DNS results with the model results of KM82 to investigate the validity 
of the model in predicting the steady states in both the regimes and also the location of bifurcation that 
leads to the transition between the two regimes. All the comparisons shown in this subsection correspond to 
e = 5 X 10“^, yd = 1 and n = 1 with a channel length = 27r. Fig. 10(a) [ shows the mean streamwise velocity 
profiles in the Hartmann regime compared to the prediction of KM82 at various values of the parameter 
Q. It can be observed that the model is accurate in this regime in terms of the magnitude of axial velocity 
although some differences can be seen in the shape of the profiles. 

However this is in contrast to the behaviour in the Poiseuille regime (see Fig. |lQ(b)| ), where the model 
strongly underpredicts the axial velocity Vx and significant differences are observed in the shape of the 
velocity profiles. These differences can be attributed to the effect of non-linearity which is more pronounced 
in the Poiseuille regime and the fact that the non-linear terms are neglected in the model equations of 
KM82. It must be pointed out that in the case of Poiseuille regime, as the final state of the flow is strongly 
time-dependent, the axial velocity profiles from the DNS are obtained by time and streamwise averaging. 

The dependance of the core axial velocity (Uc = Vx z = C) on Q is shown in Fig. ll(a)[ The 
bifurcation from the Hartmann regime to the Poiseuille regime is observed at Q ~ 0.43, which is very close 
to that predicted by KM82 and the shape of the curve is in close match. The fact that non-linearity leads 
to the differences is confirmed through a simulation that we performed dropping out the non-linear term 
(^(v • V)t;) in (|6|). Fig. |ll(a) shows that the curve obtained from DNS without the non-linear term tends 
very close to that of KM82. 
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(a) (b) 


Fig. 11. (a) Comparison of DNS results with the model predictions for e = 5 x 10“^, /^ = 1 and /^ = 1. The 
dotted line (denoted as DNS-NLT) indicates results obtained from DNS by excluding the non-linear term., 
(b) Streamwise core velocities Uc as a function of Q obtained from DNS for various values of e. Dotted lines 
indicate that no hysteresis is observed. Parameters are /d = 1 and k. = 1. 


4-.4- Effect of parameters e, /3 and n 

In this subsection, we present the effect of the parameters e, /3 and n on the nature and location of the 
bifurcation along with the magnitude of core velocity in the final state. Fig. |ll(b)| shows Uc versus Q at 
different values of e. It is clear that at higher values of e or lower Reynolds numbers, the transition from the 
Hartmann to Poiseuille regimes does not show a distict jump but rather occurs in a continuous manner. In 
the parameter space with a clear bifurcation, the value of Q at which the jump occurs is almost independent 
of the hydrodynamic Reynolds number (or e). In addition, when e is low, a two-valued solution or hysteresis 
is observed near the bifurcation point (e.g. near Q ^ 0.43 for e = 5 x 10“^), depending on the direction in 
which the steady state is approached. Such a hysteresis effect was also predicted by KM82. 

The effect of e on the core velocity {Uc) is negligible in the Hartmann regime but is strong in the Poiseuille 
zone. All these observations are akin (qualitatively) to the predictions of KM82. Interestingly, the effect of 
the magnetic Reynolds number (3 on the Uc-Q curve is very similar to that of the hydrodynamic Reynolds 
number, with higher levels of flux expulsion and core flow occuring at higher values of /3, as can be seen in 
Fig. I2(a)[ It is interesting to note that steady state solutions of KM82 are independent of the magnetic 
Reynolds number due to the association of /3 with only the non-linear term (which is neglected in the model) 
in the momentum balance, as can be seen from equations ([6]) and (HI. Furthermore, the location of the 
bifurcation to the Poiseuille regime is almost unaffected by variations in e or /3. In contrast, the jump is 
observed to be very sensitive to the streamwise wavenumber {k) of the imposed magnetic field bo. This is 
shown in Fig. |I2(b)[ indicating a clear increase in the value of Q at which the bifurcation occurs and also 
the magnitude of the jump when n is decreased. This is very similar to the dependance on n predicted by 
the inviscid version of KM82, i.e. by using e = 0 in equation ([14]). In specific, for n = 0.5, inviscid KM82 
predicts Qc ~ 0.55 as compared to Qc ~ 0-59 from viscous KM82 and Qc ~ 0.6 obtained from DNS with 
/3 = 1 and e = 5xl0“^. At a higher wavenumber {k = 2), KM82 in the inviscid limit predicts a jump at 
a value of Qc = 0.15 whereas the viscous KM82 shows a continuous transition between the two regimes. 
Interestingly in this case, DNS shows no clear demarcation between the Hartmann and the Poiseuille regimes, 
although a very small fall (rather than a jump) in Uc is observed when Q is increased from 0.17 to 0.18 and 
a corresponding onset of near wall recirculation zones at Q = 0.18. 
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Fig. 12. (a) Effect of variation of magnetic Reynolds number (3 in the Uc-Q plane, obtained from DNS. 
Parameters are hz = 1 and e = 5 x 10“^. (b) Effect of variation of streamwise wavenumber hz of the imposed 
magnetic field, obtained from DNS and KM82. Parameters are /3 = 1 and e = 5 x 10“^. 


5. Concluding remarks 

We presented results of direct numerical simulations of the dynamic runaway effect due to flux expulsion 
in a plane channel MHD flow. General features of the flow and magnetic fields in the two regimes - the 
Hartmann and Poiseuille regimes - were studied. Our results show that the one-dimensional model of Kamkar 
and Moffatt El is relatively accurate in the Hartmann regime and in the prediction of the location of the 
bifurcation to the Poiseuille regime (in the Uc-Q plane). However, significant differences in core velocity 
predictions are observed in the Poiseuille regime, attributed to the neglect of the non-linear terms in the 
Navier-Stokes equation. The Poiseuille regime is seen to be strongly unsteady similar to that of travelling 
waves, but does not show spatial irregularity for the parameters that we discussed here. The location of the 
bifurcation is found to be independent of the hydrodynamic and magnetic Reynolds numbers, but however 
is strongly affected by the wavenumber of the imposed magnetic field as lower streamwise wavenumber {k) 
leads to bifurcation to the Poiseuille regime at much higher values of Q and vice versa. Einally, in contrast 
to the model, the magnetic Reynolds number (/3) has a substantial effect on the Uc-Q curve, that is similar 
to the effect of the hydrodynamic Reynolds number (e“^). 

A full three dimensional DNS of this configuration to study well developed turbulence that might ensue 
in the Poiseuille regime can be a scope for future research. Also of interest would be to quantify the rate 
of reconnection of magnetic field lines leading to fiux expulsion in the core. This would also become much 
more demanding in three dimensions since the field line topologies are more complex. Eurthermore, in the 
case of electrically insulating walls which are of practical interest, magnetic boundary conditions could be 
non-trivial, as the secondary magnetic field transcends the channel walls and would require matching the 
exterior and interior magnetic fields at the wall boundaries. 
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